/*
    Copyright (C) 2009 Giacomo Spigler

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/


#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#include <cv.h>

#include <cutil.h>

#include "standard.h"
#include "gabor.h"


texture<unsigned char, 2, cudaReadModeElementType> c1_0;
texture<unsigned char, 2, cudaReadModeElementType> c1_1;
texture<unsigned char, 2, cudaReadModeElementType> c1_2;
texture<unsigned char, 2, cudaReadModeElementType> c1_3;





__global__ void S2_GPU(int *w, int *h, size_t pitch, int nums2, int n4, unsigned char *patch, unsigned char *s2, int psize, int basenum) {
  unsigned int scale, i;
  int x, y;
  int j, k;

  //S2 sizes
//  int w_0=(w[0]-2*GABOR_W-C1_SIZE)/C1_STEP + 1  -4; //-4 for all patches
  int h_0=(h[0]-2*GABOR_W-C1_SIZE)/C1_STEP + 1  -4;

  //C1 sizes
  int hh0=(h[0]-2*GABOR_W-C1_SIZE)/C1_STEP + 1;

  float alpha=(psize/4);
  alpha=(alpha*alpha)*2.0;


  //for scale
  for(scale=0; scale<NUMSCALES-1; scale++) {
    int w_out=(w[scale]-2*GABOR_W-C1_SIZE)/C1_STEP + 1  -psize; //patch size
    int h_out=(h[scale]-2*GABOR_W-C1_SIZE)/C1_STEP + 1  -psize;

    int sz=w_out*h_out;


    int s2coord=pitch*h_0*scale;  //h_0*scale*w_0*nums2;

    int baseY=hh0*scale;


    if(w_out>=0 && h_out>=0) {

      //for pixel
      for(i=blockIdx.x; i<sz; i+=gridDim.x) {
        x=i%w_out;
        y=i/w_out;


        //pre-load input patch in shared memory?


        int sca=psize*psize*4;
        //for patch
        for(int p=threadIdx.x; p<n4; p+=blockDim.x) {
          float val=0.0;
          int p64=p*sca;

          for(j=0; j<psize; j++) {
            int j16=j*psize*4;

            for(k=0; k<psize; k++) {
              int k4=k*4;

              float ci=tex2D(c1_0, x+j, baseY+y+k);
              float pt=patch[p64+j16+k4];
              val+=(ci-pt)*(ci-pt);

              ci=tex2D(c1_1, x+j, baseY+y+k);
              pt=patch[p64+j16+k4+1];
              val+=(ci-pt)*(ci-pt);

              ci=tex2D(c1_2, x+j, baseY+y+k);
              pt=patch[p64+j16+k4+2];
              val+=(ci-pt)*(ci-pt);

              ci=tex2D(c1_3, x+j, baseY+y+k);
              pt=patch[p64+j16+k4+3];
              val+=(ci-pt)*(ci-pt);
            }
          }



          //write output
          val/=4096.0;//4096.0;
          val/=alpha;



          //            (y*w_0+x)*nums2
          s2[s2coord + y*pitch+x*nums2 + basenum + p] = (unsigned char)( __expf(-val)*255.0);




        }


      }

    }


  }

}







extern "C" {


void ComputeS2_3(unsigned char *c1_dev[4], int w_host[NUMSCALES], int h_host[NUMSCALES], int *w_dev, int *h_dev, size_t c1_pitch, size_t s2_pitch, unsigned char *s2_dev, int nums2, int n4, int n8, unsigned char *patch4x4, unsigned char *patch8x8,   IplImage *pyr[NUMSCALES]) {
  int w_0=(w_host[0]-2*GABOR_W-C1_SIZE)/C1_STEP + 1;
  int h_0=(h_host[0]-2*GABOR_W-C1_SIZE)/C1_STEP + 1;

  cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<unsigned char>(); //(8, 0, 0, 0, 

  cudaBindTexture2D(0, &c1_0, c1_dev[0], &channelDesc, w_0, h_0*(NUMSCALES-1), c1_pitch);
  cudaBindTexture2D(0, &c1_1, c1_dev[1], &channelDesc, w_0, h_0*(NUMSCALES-1), c1_pitch);
  cudaBindTexture2D(0, &c1_2, c1_dev[2], &channelDesc, w_0, h_0*(NUMSCALES-1), c1_pitch);
  cudaBindTexture2D(0, &c1_3, c1_dev[3], &channelDesc, w_0, h_0*(NUMSCALES-1), c1_pitch);


  cudaMemset(s2_dev, 0, s2_pitch*(h_0-4)*(NUMSCALES-1) ); //


  S2_GPU<<<S2BLOCKS, S2THREADS>>>(w_dev, h_dev, s2_pitch, nums2, n4, patch4x4, s2_dev, 4, 0);
  S2_GPU<<<S2BLOCKS, S2THREADS>>>(w_dev, h_dev, s2_pitch, nums2, n8, patch8x8, s2_dev, 8, n4);
  cudaThreadSynchronize();




//TO CHECK!
/*
int scal=0;
int feat=0;

w_0-=4;
h_0-=4;
unsigned char *buf=(unsigned char *)malloc(w_0*h_0*nums2);

cvZero(pyr[0]);
cudaMemcpy2D(buf, w_0*nums2, s2_dev+scal*w_0*h_0*nums2, s2_pitch, w_0*nums2, h_0, cudaMemcpyDeviceToHost);
//for(int i=0; i<10000; i+=nums2) printf("%d\n", buf[i]);

int mmax=0;
for(int x=0; x<w_0; x++) {
  for(int y=0; y<h_0; y++) {
    pyr[0]->imageData[y*pyr[0]->width+x] = buf[(y*w_0+x)*nums2+feat];
    mmax=max(mmax, buf[(y*w_0+x)*nums2+feat]);
  }
}
printf("* max: %d\n", mmax);

*/





  cudaUnbindTexture(c1_0);
  cudaUnbindTexture(c1_1);
  cudaUnbindTexture(c1_2);
  cudaUnbindTexture(c1_3);

}









}


